Treg in JIA dataset from Adam Croft's data

suppressPackageStartupMessages({
    library(Seurat)
    library(dplyr)
    library(data.table)
    library(ggplot2)
    library(patchwork)
    library(harmony)
    library(symphony)
    library(presto)
    library(singlecellmethods)
    library(purrr)
    library(pheatmap)
    library(ggthemes)
})

set.seed(1)
source("utils.r")
source("myfun.r")

data.path <- "/data/brennerlab/Shani/projects/Treg/data/external_datasets/JIA_Croft/"

list.files(data.path)

1. 'MAPJAG.RData'
2. 'object'
3. 'umap_embeddings.rds'

metadata <- read.table(paste0(data.path, "object/meta.tsv"))
head(metadata)   
# meta.donor <- read.csv(paste0(data.path, "all.meta.CD.csv")) %>% group_by(biosample_id) %>% distinct() 
# head(meta.donor)  
dim(metadata)

A data.frame: 6 × 40

	orig.ident	nCount_RNA	nFeature_RNA	percent.mt	nCount_Protein	nFeature_Protein	T_clonotype_id	T_cdr3s_aa	B_clonotype_id	B_cdr3s_aa	⋯	global1	onset	maxjoint0623	label15	simple15	clusters2312	oldclusters	label15s	subclusters	main_clusters
	<chr>	<int>	<int>	<dbl>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	⋯	<chr>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
AAACCTGAGACAGGCT-1_1	PBMC1	2437	1009	2.379975	263	93	clonotype3221	TRB:CAISQTSATYEQYF	NA	NA	⋯	T cell	518	1	CD4+ T	T cells	CD4+ naive/central memory T	CD4+ naive	CD4+ T cells	CD4+ naive/central memory T	CD4+ T cells
AAACCTGAGACATAAC-1_1	PBMC1	3180	1428	7.106918	351	103	clonotype7162	TRB:CASSVDLAGATDTQYF	NA	NA	⋯	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells
AAACCTGAGAGGTAGA-1_1	PBMC1	3189	1486	2.759486	289	90	NA	NA	NA	NA	⋯	NK/yd T	518	1	NK cells	NK cells	IFNG+ NK	IFNG+ NK	NK cells	IFNG+ NK	NK cells
AAACCTGAGCCCTAAT-1_1	PBMC1	2057	884	4.472533	297	89	clonotype5326	TRB:CASSYLAGGTNEQFF	NA	NA	⋯	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells
AAACCTGAGCCGCCTA-1_1	PBMC1	2645	981	2.684310	292	95	clonotype6751	TRB:CASSPDRVLSAEAFF	NA	NA	⋯	T cell	518	1	CD4+ T	T cells	CD4+ naive/central memory T	CD4+ naive	CD4+ T cells	CD4+ naive/central memory T	CD4+ T cells
AAACCTGAGCGTTCCG-1_1	PBMC1	2306	951	2.385082	223	86	clonotype9493	TRB:CASSYPRERHSGANVLTF;TRA:CAVNRDDKIIF	NA	NA	⋯	T cell	518	1	Naive CD8+ T	T cells	CD8+ naive T	CD8+ naive	CD8+ T cells	CD8+ naive T	CD8+ T cells

1. 242190
2. 40

unique(metadata$main_clusters)
unique(metadata$subclusters)
unique(metadata$label15)

metadata %>% filter(subclusters == "CD4+ FOXP3+ Tregs") %>% dim
metadata %>% filter(main_clusters == "CD4+ T cells") %>% dim

Load data & organize obj

counts <- Read10X(paste0(data.path, "/object/"))

obj <- CreateSeuratObject(counts = counts,  project = "JIA_Croft", min.cells=0, min.features=0)
obj

Edit metadata

obj@meta.data%>% head

A data.frame: 6 × 3

	orig.ident	nCount_RNA	nFeature_RNA
	<fct>	<dbl>	<int>
AAACCTGAGACAGGCT-1_1	JIA_Croft	2437	1009
AAACCTGAGACATAAC-1_1	JIA_Croft	3180	1428
AAACCTGAGAGGTAGA-1_1	JIA_Croft	3189	1486
AAACCTGAGCCCTAAT-1_1	JIA_Croft	2057	884
AAACCTGAGCCGCCTA-1_1	JIA_Croft	2645	981
AAACCTGAGCGTTCCG-1_1	JIA_Croft	2306	951

mt <- obj@meta.data %>% tibble::rownames_to_column("cellID") %>% 
left_join(metadata %>% tibble::rownames_to_column("cellID"), by= join_by(cellID == cellID))

obj <- AddMetaData(obj, mt)
obj@meta.data %>% head

A data.frame: 6 × 47

	orig.ident	nCount_RNA	nFeature_RNA	cellID	orig.ident.x	nCount_RNA.x	nFeature_RNA.x	orig.ident.y	nCount_RNA.y	nFeature_RNA.y	⋯	global1	onset	maxjoint0623	label15	simple15	clusters2312	oldclusters	label15s	subclusters	main_clusters
	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>	<int>	<chr>	<int>	<int>	⋯	<chr>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
AAACCTGAGACAGGCT-1_1	JIA_Croft	2437	1009	AAACCTGAGACAGGCT-1_1	JIA_Croft	2437	1009	PBMC1	2437	1009	⋯	T cell	518	1	CD4+ T	T cells	CD4+ naive/central memory T	CD4+ naive	CD4+ T cells	CD4+ naive/central memory T	CD4+ T cells
AAACCTGAGACATAAC-1_1	JIA_Croft	3180	1428	AAACCTGAGACATAAC-1_1	JIA_Croft	3180	1428	PBMC1	3180	1428	⋯	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells
AAACCTGAGAGGTAGA-1_1	JIA_Croft	3189	1486	AAACCTGAGAGGTAGA-1_1	JIA_Croft	3189	1486	PBMC1	3189	1486	⋯	NK/yd T	518	1	NK cells	NK cells	IFNG+ NK	IFNG+ NK	NK cells	IFNG+ NK	NK cells
AAACCTGAGCCCTAAT-1_1	JIA_Croft	2057	884	AAACCTGAGCCCTAAT-1_1	JIA_Croft	2057	884	PBMC1	2057	884	⋯	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells
AAACCTGAGCCGCCTA-1_1	JIA_Croft	2645	981	AAACCTGAGCCGCCTA-1_1	JIA_Croft	2645	981	PBMC1	2645	981	⋯	T cell	518	1	CD4+ T	T cells	CD4+ naive/central memory T	CD4+ naive	CD4+ T cells	CD4+ naive/central memory T	CD4+ T cells
AAACCTGAGCGTTCCG-1_1	JIA_Croft	2306	951	AAACCTGAGCGTTCCG-1_1	JIA_Croft	2306	951	PBMC1	2306	951	⋯	T cell	518	1	Naive CD8+ T	T cells	CD8+ naive T	CD8+ naive	CD8+ T cells	CD8+ naive T	CD8+ T cells

unique(colnames(obj@meta.data))

1. 'orig.ident'
2. 'nCount_RNA'
3. 'nFeature_RNA'
4. 'cellID'
5. 'orig.ident.x'
6. 'nCount_RNA.x'
7. 'nFeature_RNA.x'
8. 'orig.ident.y'
9. 'nCount_RNA.y'
10. 'nFeature_RNA.y'
11. 'percent.mt'
12. 'nCount_Protein'
13. 'nFeature_Protein'
14. 'T_clonotype_id'
15. 'T_cdr3s_aa'
16. 'B_clonotype_id'
17. 'B_cdr3s_aa'
18. 'scrublet_score'
19. 'integrated_snn_res.0.1'
20. 'integrated_snn_res.0.15'
21. 'seurat_clusters'
22. 'integrated_snn_res.0.12'
23. 'percent.ribro'
24. 'percent.Hb'
25. 'S.Score'
26. 'G2M.Score'
27. 'Phase'
28. 'old.ident'
29. 'patient'
30. 'TYPE'
31. 'sample'
32. 'sex'
33. 'batch'
34. 'age'
35. 'condition'
36. 'njoint'
37. 'ana'
38. 'global1'
39. 'onset'
40. 'maxjoint0623'
41. 'label15'
42. 'simple15'
43. 'clusters2312'
44. 'oldclusters'
45. 'label15s'
46. 'subclusters'
47. 'main_clusters'
48. 'percent.mito'

unique(obj$global1)
unique(obj$Phase)
unique(obj$TYPE)
unique(obj$sample)
unique(obj$global1)

'T cell'

1. 'G1'
2. 'S'
3. 'G2M'

1. 'Blood'
2. 'SF'
3. 'Tissue'

1. 'B91'
2. 'B84'
3. 'B95'
4. 'B811'
5. 'B88'
6. 'B98'
7. 'B97'
8. 'B814'
9. 'B817'
10. 'SF91'
11. 'SF84'
12. 'SF95'
13. 'SF811'
14. 'SF88'
15. 'SF98'
16. 'SF97'
17. 'SF814'
18. 'SF817'
19. 'T95'
20. 'T84'
21. 'T811'
22. 'T88'
23. 'T98'
24. 'T814'
25. 'T817'
26. 'T910'
27. 'T96'
28. 'T99'

'T cell'

donor_info <- metadata %>% select(patient, TYPE, sample, sex, batch, age, condition, njoint, ana) %>% unique()
donor_info

A data.frame: 28 × 9

	patient	TYPE	sample	sex	batch	age	condition	njoint	ana
	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<int>	<chr>
AAACCTGAGACAGGCT-1_1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N
AAACCTGAGACTCGGA-1_2	84	Blood	B84	F	batch1	15.8	Poly-neg	24	N
AAACCTGAGACCACGA-1_3	95	Blood	B95	M	batch1	3.7	P-Oligo	4	N
AAACCTGAGGTGATTA-1_4	811	Blood	B811	M	batch2	10.8	P-Oligo	1	N
AAACCTGAGGATGCGT-1_5	88	Blood	B88	F	batch2	15.7	RF+Poly	23	N
AAACCTGAGACAGGCT-1_6	98	Blood	B98	F	batch2	9.3	P-Oligo	3	Y
AAACCTGAGGACCACA-1_7	97	Blood	B97	F	batch2	1.2	Poly-neg	7	Y
AAACCTGAGAGTGAGA-1_8	814	Blood	B814	M	batch3	14.8	PsA	3	NR
AAACCTGAGAGAGCTC-1_9	817	Blood	B817	F	batch3	9.3	P-Oligo	1	Y
AAACCTGAGCGTGAAC-1_10	91	SF	SF91	F	batch1	6.9	P-Oligo	1	N
AAACCTGAGACACGAC-1_11	84	SF	SF84	F	batch1	15.8	Poly-neg	24	N
AAACCTGAGCAACGGT-1_12	95	SF	SF95	M	batch1	3.7	P-Oligo	4	N
AAACCTGAGGAATGGA-1_13	811	SF	SF811	M	batch2	10.8	P-Oligo	1	N
AAACCTGAGAGCTATA-1_14	88	SF	SF88	F	batch2	15.7	RF+Poly	23	N
AAACCTGAGCCACTAT-1_15	98	SF	SF98	F	batch2	9.3	P-Oligo	3	Y
AAACCTGAGAATCTCC-1_16	97	SF	SF97	F	batch2	1.2	Poly-neg	7	Y
AAACCTGTCAGAGACG-1_17	814	SF	SF814	M	batch3	14.8	PsA	3	NR
AAACCTGAGAAGGTTT-1_18	817	SF	SF817	F	batch3	9.3	P-Oligo	1	Y
AAACCTGAGAAGGACA-1_19	95	Tissue	T95	M	batch2	3.7	P-Oligo	4	N
AAACCTGAGAATTCCC-1_20	84	Tissue	T84	F	batch2	15.8	Poly-neg	24	N
AAACCTGAGTTTCCTT-1_21	811	Tissue	T811	M	batch2	10.8	P-Oligo	1	N
AAACCTGAGACCTAGG-1_22	88	Tissue	T88	F	batch2	15.7	RF+Poly	23	N
AAACCTGAGAAGGTGA-1_23	98	Tissue	T98	F	batch2	9.3	P-Oligo	3	Y
AAACCTGAGCAATATG-1_24	814	Tissue	T814	M	batch3	14.8	PsA	3	NR
AAACCTGAGAGAGCTC-1_25	817	Tissue	T817	F	batch3	9.3	P-Oligo	1	Y
AAACCTGAGACTCGGA-1_26	910	Tissue	T910	F	batch3	3.0	Ex-Oligo	1	Y
AAACCTGAGAATTCCC-1_27	96	Tissue	T96	F	batch3	4.8	PsA	1	N
AAACCTGAGAGGACGG-1_28	99	Tissue	T99	F	batch3	2.6	P-Oligo	2	NR

table(donor_info$batch, donor_info$patient)
table(donor_info$patient, donor_info$TYPE)

        
         84 88 91 95 96 97 98 99 811 814 817 910
  batch1  2  0  2  2  0  0  0  0   0   0   0   0
  batch2  1  3  0  1  0  2  3  0   3   0   0   0
  batch3  0  0  0  0  1  0  0  1   0   3   3   1

     
      Blood SF Tissue
  84      1  1      1
  88      1  1      1
  91      1  1      0
  95      1  1      1
  96      0  0      1
  97      1  1      0
  98      1  1      1
  99      0  0      1
  811     1  1      1
  814     1  1      1
  817     1  1      1
  910     0  0      1

# Keep T cells only
Idents(obj) <- "global1"
obj <- subset(obj, ident = "T cell")

Basic QC & preproccessing

rownames(obj)[grep(pattern = "^MT-", rownames(obj))]

1. 'MT-ND1'
2. 'MT-ND2'
3. 'MT-CO1'
4. 'MT-CO2'
5. 'MT-ATP8'
6. 'MT-ATP6'
7. 'MT-CO3'
8. 'MT-ND3'
9. 'MT-ND4L'
10. 'MT-ND4'
11. 'MT-ND5'
12. 'MT-ND6'
13. 'MT-CYB'

obj[["percent.mito"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
VlnPlot(obj, features = "percent.mito", pt.size = 0.0, group.by = "patient") + theme(legend.position = "none") + ggtitle(NULL) +
  xlab(NULL) + ylab("% mitochondrial genes") + geom_hline(yintercept=15)

VlnPlot(obj, features = "percent.mito", pt.size = 0.0, group.by = "batch") + theme(legend.position = "none") + ggtitle(NULL) +
  xlab(NULL) + ylab("% mitochondrial genes") + geom_hline(yintercept=15)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

#Gene number
VlnPlot(obj, features = "nFeature_RNA", pt.size = 0.0, group.by = "patient") + theme(legend.position = "none") + ggtitle(NULL) +
  xlab(NULL) + ylab("Number of Genes")
VlnPlot(obj, features = "nFeature_RNA", pt.size = 0.0, group.by = "batch") + theme(legend.position = "none") + ggtitle(NULL) +
  xlab(NULL) + ylab("Number of Genes")

VlnPlot(obj, features = "nCount_RNA", pt.size = 0.0, group.by = "patient") + theme(legend.position = "none") + ggtitle(NULL) +
  xlab(NULL) + ylab("nUMI")
VlnPlot(obj, features = "nCount_RNA", pt.size = 0.0, group.by = "batch") + theme(legend.position = "none") + ggtitle(NULL) +
  xlab(NULL) + ylab("nUMI")

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

fig.size(8,8)
suppressMessages(
    obj@meta.data %>% 
  	ggplot(aes(x=nCount_RNA, y=nFeature_RNA, color=percent.mito)) + 
  	geom_point() + 
	scale_colour_gradient(low = "gray90", high = "black") +
    facet_wrap(~orig.ident) +
  	stat_smooth(method=lm) +
  	scale_x_log10() + 
  	scale_y_log10() + 
  	theme_classic() +
  	geom_vline(xintercept = 500) +
  	geom_hline(yintercept = 250) +
    theme(text = element_text(size = 14)) 
    )

#Data already filtered
# obj <- subset(obj, subset = nCount_RNA > 500 & nFeature_RNA > 250 &  percent.mito < 20)

#initial dimensionality reduction on all samples to identify contaminating clusters
obj <- NormalizeData(obj, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>% 
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>% 
    ScaleData() %>% 
    RunPCA(verbose = F)

Centering and scaling data matrix

set.seed(1)
obj <- RunHarmony(obj, group.by.vars=c("batch", "patient", "TYPE"), reduction.use = "pca", 
                 plot_convergence = TRUE, max_iter = 10,
                 early_stop = F)

Transposing data matrix

Initializing state using k-means centroids initialization

Harmony 1/10

Harmony 2/10

Harmony 3/10

Harmony 4/10

Harmony 5/10

Harmony 6/10

Harmony 7/10

Harmony 8/10

Harmony 9/10

Harmony 10/10

obj

An object of class Seurat 
30803 features across 100098 samples within 1 assay 
Active assay: RNA (30803 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 2 dimensional reductions calculated: pca, harmony

set.seed(1)
obj <- Run_uwot_umap(obj)
obj <- FindClusters(obj, graph.name = 'humap_fgraph', resolution = 1, verbose = TRUE)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 100098
Number of edges: 1234828

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8233
Number of communities: 14
Elapsed time: 14 seconds

fig.size(8,8)
DimPlot(obj, label = T, repel = T, pt.size = 0.5, label.size = 6) + #scale_color_tableau("Tableau 20") +
   theme(text = element_text(size = 20))
DimPlot(obj, group.by = "batch")
DimPlot(obj, group.by = "patient")
DimPlot(obj, group.by = "TYPE")
DimPlot(obj, group.by = "condition")

fig.size(8,16)
DimPlot(obj, group.by = "subclusters")

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

fig.size(12, 20)
FeaturePlot(obj, c("CD3E", "CD4", "FOXP3", "IL2RA", "CTLA4", "IKZF2", "CD8A", "CXCR6"), ncol = 4, raster=T)
FeaturePlot(obj, c("CD3E", "CD4", "FOXP3", "IL2RA", "CTLA4", "IKZF2", "CD8A", "CXCR6"), ncol = 4, order = T, raster=T)

fig.size(4,10)
VlnPlot(obj, "CD4", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CD8A", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CD3E", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CD3D", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "FOXP3", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "IL2RA", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CXCL13", group.by = "seurat_clusters", raster=T)

table(obj$subclusters, obj$seurat_clusters)

                             
                                 0    1    2    3    4    5    6    7    8    9
  Activated NK-like T          107    7    9   10 1034   10 7632    1    1    6
  CD4+ FOXP3+ Tregs              3   38   97   12    3   37    4  113 4327    4
  CD4+ KLRB1+ memory T         303 8761  654   71   40  158   42  527   69   76
  CD4+ naive/central memory T   30 3670 9539  578  253 9681   34  112   32    8
  CD56+bright NK                 1    4    1    0    0    1    8    0    0    0
  CD8+ GZMB+/GZMK+ memory T   7870    6    1   26    1    1   62   12    3 1548
  CD8+ GZMK+ memory T         5853  125  186  105  756   12  286   21    3  215
  CD8+ naive T                 113   40  708  109 8401  278  354    4    0   13
  CXCL13+ TpH                   72 1774  333  161   11  157   18 5327 1209  449
  Cycling Myeloid                0    0    1    0    0    0    8    0    0    0
  Cycling T                      4    7    0    1    2    0   50    1    1    3
  IFNG+ NK                      72    0    0    0    0    0    5    0    0    3
  MAIT cells                     3   33    4    3    0    0    0    0    0    0
  NK-like ILC1                  12    0   14   79    1    0    1    0    0    0
  THEMIS+ IL7R+ ILC            235   73  635 9945   67  197   62   46   82   30
  Vdelta2 yd                    27    4    0    0    1    0   10    0    0    0
                             
                                10   11   12   13
  Activated NK-like T           14    0    0   19
  CD4+ FOXP3+ Tregs              1   80   13    6
  CD4+ KLRB1+ memory T           0    0  131   40
  CD4+ naive/central memory T    3    7  136   64
  CD56+bright NK                 2    0    0    9
  CD8+ GZMB+/GZMK+ memory T      4    0   10   16
  CD8+ GZMK+ memory T           24    0   53   23
  CD8+ naive T                   7    0   57   22
  CXCL13+ TpH                    2   29   90   17
  Cycling Myeloid                0    0    0    0
  Cycling T                      0    0    0    0
  IFNG+ NK                       5    0    0   34
  MAIT cells                     0    0    0    0
  NK-like ILC1                1183    0    0   66
  THEMIS+ IL7R+ ILC            367  990    4  112
  Vdelta2 yd                     0    0    0    5

subset Tregs I

#subset 
tregs <- subset(obj, ident = c(8))
tregs

An object of class Seurat 
30803 features across 5727 samples within 1 assay 
Active assay: RNA (30803 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 3 dimensional reductions calculated: pca, harmony, humap

redo preproccessing

tregs <- NormalizeData(tregs, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>% 
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>% 
    ScaleData() %>% 
    RunPCA(verbose = F)

Centering and scaling data matrix

set.seed(1)
tregs <- RunHarmony(tregs, group.by.vars=c("batch", "patient", "TYPE"), reduction.use = "pca", 
                 plot_convergence = TRUE, max_iter = 10,
                 early_stop = T)

Transposing data matrix

Initializing state using k-means centroids initialization

Harmony 1/10

Harmony 2/10

Harmony 3/10

Harmony 4/10

Harmony converged after 4 iterations

set.seed(1)
tregs <- Run_uwot_umap(tregs)
tregs <- FindClusters(tregs, graph.name = 'humap_fgraph', resolution = 0.8, verbose = TRUE)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 5727
Number of edges: 69517

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.6827
Number of communities: 9
Elapsed time: 0 seconds

fig.size(8,8)
DimPlot(tregs, label = T, repel = T, pt.size = 0.5, label.size = 6) + 
  scale_color_tableau("Tableau 20") + theme(text = element_text(size = 20))
DimPlot(tregs, group.by = "batch")
DimPlot(tregs, group.by = "patient")
DimPlot(tregs, group.by = "TYPE")
DimPlot(tregs, group.by = "condition")

fig.size(8,16)
DimPlot(tregs, group.by = "subclusters")

fig.size(12, 15)
FeaturePlot(tregs, c("CD3E", "CD4", "FOXP3", "IL2RA", "CTLA4", "IKZF2", "AREG", "TCF7", "CXCR6"), ncol = 3)

fig.size(4,10)
VlnPlot(tregs, "FOXP3", group.by = "seurat_clusters")
VlnPlot(tregs, "IL2RA", group.by = "seurat_clusters")
VlnPlot(tregs, "AREG", group.by = "seurat_clusters")
VlnPlot(tregs, "TCF7", group.by = "seurat_clusters")

# comapre to source annotations
table(tregs$subclusters, tregs$seurat_clusters)

                             
                                0   1   2   3   4   5   6   7   8
  Activated NK-like T           0   0   0   0   0   0   1   0   0
  CD4+ FOXP3+ Tregs           629 903 872 576 802 362 123  35  25
  CD4+ KLRB1+ memory T         27  18   5  12   0   1   6   0   0
  CD4+ naive/central memory T  28   3   0   1   0   0   0   0   0
  CD8+ GZMB+/GZMK+ memory T     0   0   0   0   0   0   3   0   0
  CD8+ GZMK+ memory T           0   0   0   1   0   1   1   0   0
  CXCL13+ TpH                 455 213  99 283   8 103  41   5   2
  Cycling T                     0   0   0   0   1   0   0   0   0
  THEMIS+ IL7R+ ILC            11   0   1   0   2  66   2   0   0

cluster DEGs

Tregs.markers <- wilcoxauc(tregs, "seurat_clusters")

options(repr.matrix.max.cols=150, repr.matrix.max.rows=200)

Tregs.markers %>% 
    group_by(group) %>% 
    filter(padj < 0.05 & auc > 0.6) %>% 
    arrange(group, desc(logFC)) %>% 
    #slice_max(n = 40, order_by = logFC)%>% 
    dplyr::select(feature, group) %>% 
    mutate(row = row_number()) %>% 
    tidyr::pivot_wider(names_from = group, values_from = feature) -> Treg.show
Treg.show[1:50,]

A tibble: 50 × 10

row	0	1	2	3	4	5	6	7	8
<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
1	LEF1	KLF2	LGALS1	FOS	TNFRSF18	MTRNR2L12	CCL5	NR4A1	TRBV7-6
2	RPS12	EMP3	HLA-DRB1	JUNB	TNFRSF4	NEAT1	NKG7	DUSP2	TRBV7-4
3	EEF1B2	TXNIP	S100A4	DUSP1	TIGIT	SYNE2	CD8A	CD69	TRBV7-2
4	TCF7	TLE5	CD74	CREM	LINC01943	MT-ATP6	CTSW	FOS	TRBV7-9
5	RPS3A	CD52	S100A11	ZNF331	CTLA4	MALAT1	CD8B	NFKBID	TRAV38-2DV8
6	RPS13	LIME1	LGALS3	FOSB	BATF	MT-CO1	GZMA	TNFRSF4	KLRB1
7	RPL32	SELPLG	S100A6	DUSP2	CTSC	ARHGAP15	GZMK	ZFP36L1	GZMA
8	RPS8	S100A10	VIM	ZFP36	CARD16	FOXP1	CST7	CTLA4	STX10
9	RPS5	S100A6	CRIP1	DNAJB1	IL2RA	RNF213	LAG3	ZFP36	PIK3CG
10	RPS18	SH3BGRL3	S100A10	TNFAIP3	FOXP3	MT-ND5	GNLY	PIM3	CLSTN3
11	RPL13	TMSB10	ANXA2	RGCC	CD74	MT-CO3	KLRK1	RGCC	FBXL14
12	RPL30	B2M	HLA-DRA	CXCR4	SPOCK2	SKAP1	CCL4	EGR3	PRKAB1
13	TXNIP	HLA-C	CLIC1	NR4A2	SRGN	CASK	GZMB	EGR1	NA
14	RPL22	RPL14	ANXA1	JUND	IL2RB	SMCHD1	ZFP36L2	NR4A3	NA
15	RPL3	RPL18A	UCP2	CD69	TYMP	LINC02694	ANXA1	PHLDA1	NA
16	RPS27	NA	HLA-DPA1	JUN	GAPDH	MT-ND6	CD96	BATF	NA
17	RPL39	NA	MYL6	DNAJA1	TNFRSF1B	MBNL1	THEMIS	TNF	NA
18	MAL	NA	TIMP1	LMNA	ENO1	MT-CYB	CD7	SH2D2A	NA
19	RPS6	NA	EMP3	YPEL5	DUSP4	MT-CO2	GZMM	ICOS	NA
20	CCR7	NA	CD99	SLC2A3	CXCR6	PHACTR2	PDE3B	KLF6	NA
21	RPL21	NA	GZMA	TENT5C	PKM	CHST11	AOAH	BTG2	NA
22	SELL	NA	HLA-DPB1	RGS1	LAYN	EPB41	PARP8	NFKBIA	NA
23	RPS23	NA	LY6E	PPP1R15A	SAT1	MT-ND1	HCST	SRGN	NA
24	RPS21	NA	SH3BGRL3	BTG2	CST7	LPP	PTMS	CSRNP1	NA
25	RPL34	NA	HLA-DRB5	DUSP4	TPI1	MT-ND4L	GZMH	GEM	NA
26	RPLP2	NA	GAPDH	SELENOK	PHACTR2	USP15	STK17A	BIRC3	NA
27	RPL18A	NA	CD2	KLF6	CD2	RABGAP1L	S100A6	TNFRSF18	NA
28	GIMAP7	NA	TXN	PHLDA1	ACP5	AKAP13	SP100	JUND	NA
29	RPL5	NA	CTSC	FTH1	CD7	GLCCI1	LYAR	IER2	NA
30	RPS29	NA	TYMP	TSC22D3	IGFLR1	CAMK4	FXYD2	CREM	NA
31	RPL11	NA	CKLF	HSPH1	PHLDA1	PDE3B	ADGRE5	NR4A2	NA
32	RPL9	NA	ARPC1B	SRGN	MAGEH1	DPYD	CD99	TNFRSF1B	NA
33	RPL10A	NA	COTL1	UBE2S	PRDM1	GPRIN3	CBLB	DNAJA1	NA
34	RPS25	NA	IL2RG	HSPA8	CACYBP	LDLRAD4	CD81	JUNB	NA
35	RPL35A	NA	ACTG1	NR4A1	AC017002.3	TNIK	PRF1	PTPN7	NA
36	RPS27A	NA	ENO1	HSPA1B	CD27	IQGAP2	TCF25	GNG2	NA
37	RPL29	NA	CORO1A	HSP90AA1	IL2RG	MT-ND3	PIK3R5	REL	NA
38	RPS14	NA	PLP2	SOCS1	DNPH1	PIK3R5	PITPNC1	TBC1D4	NA
39	RPL36	NA	MYL12A	BTG1	PGAM1	PACS1	RUNX3	ARID5B	NA
40	RPL37	NA	AHNAK	SAMSN1	STK17B	KLF12	LITAF	NAMPT	NA
41	RPL19	NA	TPM3	TAGAP	GBP5	LRBA	APOBEC3G	RNF19A	NA
42	KLF2	NA	HCST	H3F3B	LGALS3	PRKCH	CYBA	HSP90AA1	NA
43	RPL36A	NA	CD52	HSP90AB1	ICOS	MT-ND4	ITGAE	RILPL2	NA
44	RPS3	NA	PKM	CCNH	MYL6	MT-ND2	ID2	PPP1R15A	NA
45	RPL10	NA	ACTB	NR4A3	CBLB	SAMHD1	CTSD	NEDD9	NA
46	RPL38	NA	CXCR3	MYADM	IL12RB2	NA	CXCR6	NFKB1	NA
47	RPL8	NA	HLA-DQB1	PTPN22	RNF213	NA	TLN1	ARF6	NA
48	EEF1G	NA	CYBA	IRF1	GBP2	NA	MBP	MTRNR2L12	NA
49	RPL12	NA	LSP1	HSPE1	PTPRJ	NA	EMB	STAT3	NA
50	EEF1A1	NA	RNF213	NAMPT	LAG3	NA	GSTP1	RAB8B	NA

table(tregs$seurat_clusters)

   0    1    2    3    4    5    6    7    8    9   10   11 
1538 1224 1178 1102 1018  970  576  506  434  233  152  130 

Subset Tregs II

consistent with the strategy for other datasets

# 6 is CD8
# 8 is contamination TRBV genes
tregs <- subset(tregs, ident = c(6,8), invert = T)

Final Tregs - Normalization, DimRed, Clustering

#initial dimensionality reduction on all samples to identify contaminating clusters
tregs <- NormalizeData(tregs, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>% 
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>% 
    ScaleData() %>% 
    RunPCA(verbose = F)

Centering and scaling data matrix

tregs

An object of class Seurat 
30803 features across 5523 samples within 1 assay 
Active assay: RNA (30803 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 3 dimensional reductions calculated: pca, harmony, humap

tregs@meta.data %>% head

A data.frame: 6 × 50

	orig.ident	nCount_RNA	nFeature_RNA	cellID	orig.ident.x	nCount_RNA.x	nFeature_RNA.x	orig.ident.y	nCount_RNA.y	nFeature_RNA.y	percent.mt	nCount_Protein	nFeature_Protein	T_clonotype_id	T_cdr3s_aa	B_clonotype_id	B_cdr3s_aa	scrublet_score	integrated_snn_res.0.1	integrated_snn_res.0.15	seurat_clusters	integrated_snn_res.0.12	percent.ribro	percent.Hb	S.Score	G2M.Score	Phase	old.ident	patient	TYPE	sample	sex	batch	age	condition	njoint	ana	global1	onset	maxjoint0623	label15	simple15	clusters2312	oldclusters	label15s	subclusters	main_clusters	percent.mito	humap_fgraph_res.1	humap_fgraph_res.0.8
	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>	<int>	<chr>	<int>	<int>	<dbl>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<int>	<int>	<fct>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<chr>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<int>	<chr>	<chr>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<fct>	<fct>
AAAGATGAGCAGCGTA-1_1	JIA_Croft	1350	760	AAAGATGAGCAGCGTA-1_1	JIA_Croft	1350	760	PBMC1	1350	760	3.1111111	340	99	clonotype4891	TRB:CASSLAIRNYGYTF	NA	NA	0.06086957	1	1	0	1	0.1507320	0	0.01850598	0.068792887	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells	3.1111111	8	0
AAAGATGTCCTAGTGA-1_1	JIA_Croft	1485	773	AAAGATGTCCTAGTGA-1_1	JIA_Croft	1485	773	PBMC1	1485	773	1.4814815	292	93	clonotype77	TRA:CAVRLSWGKLQF	NA	NA	0.08163265	1	1	1	1	0.1495413	0	-0.04059452	0.008258175	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.4814815	8	1
AAAGTAGCACTAGTAC-1_1	JIA_Croft	1567	847	AAAGTAGCACTAGTAC-1_1	JIA_Croft	1567	847	PBMC1	1567	847	2.1059349	330	100	clonotype1715	TRB:CASSSQNRAEQYF	NA	NA	0.05263158	1	1	0	1	0.1298440	0	0.01361041	0.082733395	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	2.1059349	8	0
AACCATGGTCGCATCG-1_1	JIA_Croft	2589	1048	AACCATGGTCGCATCG-1_1	JIA_Croft	2589	1048	PBMC1	2589	1048	0.9269988	278	89	clonotype751	TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF	NA	NA	0.15625000	1	1	0	1	0.1489572	0	0.04661555	-0.045203442	S	S	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	0.9269988	8	0
AACCATGGTTGCTCCT-1_1	JIA_Croft	2049	1102	AACCATGGTTGCTCCT-1_1	JIA_Croft	2049	1102	PBMC1	2049	1102	1.2201074	279	90	clonotype9315	TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF	NA	NA	0.03914591	1	1	1	1	0.1104492	0	-0.04041617	-0.052987582	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.2201074	8	1
AACGTTGAGGATGGAA-1_1	JIA_Croft	4757	1769	AACGTTGAGGATGGAA-1_1	JIA_Croft	4757	1769	PBMC1	4757	1769	1.8499054	293	92	clonotype1286	TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF	NA	NA	0.09497207	1	1	0	1	0.1185773	0	-0.02773367	-0.003626403	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	1.8499054	8	0

fig.size(4,4)
set.seed(1)
tregs <- RunHarmony(tregs, group.by.vars=c("batch", "patient", "TYPE"), reduction.use = "pca", 
                 plot_convergence = TRUE, max_iter = 10,
                 early_stop = T)

# ElbowPlot(obj, ndims = 50, reduction = "harmony") 
# ElbowPlot(obj, ndims = 50, reduction = "pca")

Transposing data matrix

Initializing state using k-means centroids initialization

Harmony 1/10

Harmony 2/10

Harmony 3/10

Harmony 4/10

Harmony converged after 4 iterations

set.seed(1)
tregs <- Run_uwot_umap(tregs)
tregs <- FindClusters(tregs, graph.name = 'humap_fgraph', resolution = 0.8, verbose = TRUE)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 5523
Number of edges: 67522

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.6758
Number of communities: 8
Elapsed time: 0 seconds

fig.size(6,6)
DimPlot(tregs, label = T, repel = T, pt.size = 0.5, label.size = 6) + 
  scale_color_tableau("Tableau 20") + theme(text = element_text(size = 20))
DimPlot(tregs, group.by = "batch")
DimPlot(tregs, group.by = "patient")
DimPlot(tregs, group.by = "TYPE")
DimPlot(tregs, group.by = "condition")

fig.size(8,16)
DimPlot(tregs, group.by = "subclusters")

fig.size(12, 12)
FeaturePlot(tregs, c("FOXP3", "IL2RA", "CTLA4", "IKZF2", "CXCR6", "AREG", "MKI67", "TCF7", "CXCL13"), ncol = 3)
FeaturePlot(tregs, order = T, c("FOXP3", "IL2RA", "CTLA4", "IKZF2", "CXCR6", "AREG", "MKI67", "TCF7", "CXCL13"), ncol = 3)

fig.size(5, 30)
FeaturePlot(tregs, order = T, c("FOXP3"), split.by = "seurat_clusters", ncol = 4)
FeaturePlot(tregs, order = F, c("FOXP3"), split.by = "seurat_clusters", ncol = 4)

fig.size(2,10)
VlnPlot(tregs, "FOXP3", group.by = "seurat_clusters")
VlnPlot(tregs, "IL2RA", group.by = "seurat_clusters")
VlnPlot(tregs, "CD4", group.by = "seurat_clusters")
VlnPlot(tregs, "CD8A", group.by = "seurat_clusters")
VlnPlot(tregs, "CXCL13", group.by = "seurat_clusters")

Save Obj.

saveRDS(tregs, "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs.RDS")
# treg <- readRDS("/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/CD_Tregs.RDS")

Map to RA tregs with symphony

ref <- "/data/brennerlab/Shani/projects/Treg/analysis/integrated/integrated_Tregs_SymphonyRef.rds"
tregRef<- readRDS(ref)

query = mapQuery(tregs@assays$RNA@counts,             # query gene expression (genes x cells)
                 tregs@meta.data,        # query metadata (cells x attributes)
                 tregRef,             # Symphony reference object 
                 vars = c("batch", "patient", "TYPE"),           #"PubID", "Site", "Type" # Query batch variables to harmonize over (NULL treats query as one batch)
                 do_normalize = TRUE,  # perform log(CP10k) normalization on query (set to FALSE if already normalized)
                 do_umap = TRUE)        # project query cells into reference UMAP

Normalizing

Scaling and synchronizing query gene expression

Found 1940 out of 2000 reference variable genes in query dataset

Project query cells using reference gene loadings

Clustering query cells to reference centroids

Correcting query batch effects

UMAP

All done!

summary(query)

          Length    Class      Mode   
exp       170124969 dgCMatrix  S4     
meta_data        54 data.frame list   
Z            276150 -none-     numeric
Zq_pca       276150 -none-     numeric
R            552300 -none-     numeric
Xq            44184 -none-     numeric
umap          11046 -none-     numeric

fig.size(4,6)
qp <- query$umap %>% as.data.frame %>% cbind(query$meta_data) %>% 
ggplot(aes(x=UMAP1, y=UMAP2, color=seurat_clusters)) + geom_point(size = 0.2) + 
theme_bw(base_size = 18) 

qp
fig.size(4,10)
qp + facet_wrap(facets = vars(TYPE))
qp + facet_wrap(facets = vars(patient))

# Predict query cell types
query <- knnPredict(query, tregRef, 
                   train_labels = tregRef$meta_data$cell.states, 
                   k = 30, 
                   confidence = TRUE) # calculate prediction confidence

query$meta_data %>% head

A data.frame: 6 × 56

	orig.ident	nCount_RNA	nFeature_RNA	cellID	orig.ident.x	nCount_RNA.x	nFeature_RNA.x	orig.ident.y	nCount_RNA.y	nFeature_RNA.y	percent.mt	nCount_Protein	nFeature_Protein	T_clonotype_id	T_cdr3s_aa	B_clonotype_id	B_cdr3s_aa	scrublet_score	integrated_snn_res.0.1	integrated_snn_res.0.15	seurat_clusters	integrated_snn_res.0.12	percent.ribro	percent.Hb	S.Score	G2M.Score	Phase	old.ident	patient	TYPE	sample	sex	batch	age	condition	njoint	ana	global1	onset	maxjoint0623	label15	simple15	clusters2312	oldclusters	label15s	subclusters	main_clusters	percent.mito	humap_fgraph_res.1	humap_fgraph_res.0.8	cell.states	cell.states.knn.prob	donorID	tissue	cell_type_pred_knn	cell_type_pred_knn_prob
	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>	<int>	<chr>	<int>	<int>	<dbl>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<int>	<int>	<fct>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<chr>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<int>	<chr>	<chr>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<fct>	<fct>	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>
AAAGATGAGCAGCGTA-1_1	JIA_Croft	1350	760	AAAGATGAGCAGCGTA-1_1	JIA_Croft	1350	760	PBMC1	1350	760	3.1111111	340	99	clonotype4891	TRB:CASSLAIRNYGYTF	NA	NA	0.06086957	1	1	0	1	0.1507320	0	0.01850598	0.068792887	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells	3.1111111	8	0	Naive Treg	0.4333333	91	Blood	Naive Treg	0.4333333
AAAGATGTCCTAGTGA-1_1	JIA_Croft	1485	773	AAAGATGTCCTAGTGA-1_1	JIA_Croft	1485	773	PBMC1	1485	773	1.4814815	292	93	clonotype77	TRA:CAVRLSWGKLQF	NA	NA	0.08163265	1	1	0	1	0.1495413	0	-0.04059452	0.008258175	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.4814815	8	0	Naive Treg	0.5666667	91	Blood	CD25int Treg	0.5666667
AAAGTAGCACTAGTAC-1_1	JIA_Croft	1567	847	AAAGTAGCACTAGTAC-1_1	JIA_Croft	1567	847	PBMC1	1567	847	2.1059349	330	100	clonotype1715	TRB:CASSSQNRAEQYF	NA	NA	0.05263158	1	1	1	1	0.1298440	0	0.01361041	0.082733395	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	2.1059349	8	1	Naive Treg	0.6000000	91	Blood	Naive Treg	0.6333333
AACCATGGTCGCATCG-1_1	JIA_Croft	2589	1048	AACCATGGTCGCATCG-1_1	JIA_Croft	2589	1048	PBMC1	2589	1048	0.9269988	278	89	clonotype751	TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF	NA	NA	0.15625000	1	1	1	1	0.1489572	0	0.04661555	-0.045203442	S	S	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	0.9269988	8	1	Naive Treg	0.9333333	91	Blood	Naive Treg	0.9333333
AACCATGGTTGCTCCT-1_1	JIA_Croft	2049	1102	AACCATGGTTGCTCCT-1_1	JIA_Croft	2049	1102	PBMC1	2049	1102	1.2201074	279	90	clonotype9315	TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF	NA	NA	0.03914591	1	1	6	1	0.1104492	0	-0.04041617	-0.052987582	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.2201074	8	6	CD25int Treg	0.5666667	91	Blood	ISG Treg	0.5333333
AACGTTGAGGATGGAA-1_1	JIA_Croft	4757	1769	AACGTTGAGGATGGAA-1_1	JIA_Croft	4757	1769	PBMC1	4757	1769	1.8499054	293	92	clonotype1286	TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF	NA	NA	0.09497207	1	1	1	1	0.1185773	0	-0.02773367	-0.003626403	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	1.8499054	8	1	Naive Treg	0.7666667	91	Blood	Naive Treg	0.7333333

fig.size(8, 12)
qp <- query$umap %>% as.data.frame %>% cbind(query$meta_data) %>% 
ggplot(aes(x=UMAP1, y=UMAP2, color=cell_type_pred_knn)) + geom_point(size = 0.8) + 
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) + 
scale_color_manual(values = cell.state.colors$cell.states) + 
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +coord_fixed()

qp
fig.size(10,18)
qp + facet_wrap(facets = vars(cell_type_pred_knn))
fig.size(6, 16)
qp + facet_wrap(facets = vars(TYPE))
qp + facet_wrap(facets = vars(batch))
qp + facet_wrap(facets = vars(patient))
qp + facet_wrap(facets = vars(age))

query$meta_data %>% select(subclusters, cell_type_pred_knn, cell_type_pred_knn_prob) %>% arrange(cell_type_pred_knn_prob) %>% head
table(query$meta_data$subclusters, query$meta_data$cell_type_pred_knn)

A data.frame: 6 × 3

	subclusters	cell_type_pred_knn	cell_type_pred_knn_prob
	<chr>	<fct>	<dbl>
GTGCTTCCATCGGACC-1_12	CD4+ KLRB1+ memory T	AREG+ Treg	0.2333333
ACATACGCAAGCTGAG-1_14	CXCL13+ TpH	CD25int Treg	0.2333333
TTCGAAGAGCGAAGGG-1_14	CD4+ FOXP3+ Tregs	CD25highCXCR6+ Treg	0.2580645
TCAGCTCTCGTCTGAA-1_10	CXCL13+ TpH	CD161+mem. Treg	0.2666667
CCAGCGATCTGGCGAC-1_11	CXCL13+ TpH	CD25high Treg	0.2666667
CATATGGCAGTTTACG-1_14	CD4+ FOXP3+ Tregs	AREG+ Treg	0.2666667

                             
                              Naive Treg CD25int Treg CD25high Treg
  CD4+ FOXP3+ Tregs                  285          516          1557
  CD4+ KLRB1+ memory T                 8           18            11
  CD4+ naive/central memory T         15            8             0
  CD8+ GZMK+ memory T                  0            0             0
  CXCL13+ TpH                         98          271           226
  Cycling T                            0            0             0
  THEMIS+ IL7R+ ILC                   16           11            17
                             
                              CD25highCXCR6+ Treg AREG+ Treg TNFAIP3+ Treg
  CD4+ FOXP3+ Tregs                           990        212           358
  CD4+ KLRB1+ memory T                          1         15             2
  CD4+ naive/central memory T                   0          8             0
  CD8+ GZMK+ memory T                           0          1             1
  CXCL13+ TpH                                  59        332            94
  Cycling T                                     0          0             0
  THEMIS+ IL7R+ ILC                            20         11             0
                             
                              CD161+mem. Treg ISG Treg GZM Treg Prolif.
  CD4+ FOXP3+ Tregs                        61      179       14       7
  CD4+ KLRB1+ memory T                      5        0        3       0
  CD4+ naive/central memory T               1        0        0       0
  CD8+ GZMK+ memory T                       0        0        0       0
  CXCL13+ TpH                              41       38        6       1
  Cycling T                                 0        0        0       1
  THEMIS+ IL7R+ ILC                         2        2        0       1

fig.size(8, 12)
qp + facet_wrap(facets = vars(subclusters, TYPE))
# qp + facet_wrap(facets = vars(orig.ident, Layer))

mt<- tregs@meta.data %>% tibble::rownames_to_column("rownames") %>% 
left_join(query$meta_data %>% tibble::rownames_to_column("rownames"), by = join_by(rownames))
mt %>% head

A data.frame: 6 × 111

	rownames	orig.ident.x.x	nCount_RNA.x.x	nFeature_RNA.x.x	cellID.x	orig.ident.x.x.x	nCount_RNA.x.x.x	nFeature_RNA.x.x.x	orig.ident.y.x	nCount_RNA.y.x	nFeature_RNA.y.x	percent.mt.x	nCount_Protein.x	nFeature_Protein.x	T_clonotype_id.x	T_cdr3s_aa.x	B_clonotype_id.x	B_cdr3s_aa.x	scrublet_score.x	integrated_snn_res.0.1.x	integrated_snn_res.0.15.x	seurat_clusters.x	integrated_snn_res.0.12.x	percent.ribro.x	percent.Hb.x	S.Score.x	G2M.Score.x	Phase.x	old.ident.x	patient.x	TYPE.x	sample.x	sex.x	batch.x	age.x	condition.x	njoint.x	ana.x	global1.x	onset.x	maxjoint0623.x	label15.x	simple15.x	clusters2312.x	oldclusters.x	label15s.x	subclusters.x	main_clusters.x	percent.mito.x	humap_fgraph_res.1.x	humap_fgraph_res.0.8.x	cell.states.x	cell.states.knn.prob.x	donorID.x	tissue.x	orig.ident.y.y	nCount_RNA.y.y	nFeature_RNA.y.y	cellID.y	orig.ident.x.y	nCount_RNA.x.y	nFeature_RNA.x.y	orig.ident.y.y.y	nCount_RNA.y.y.y	nFeature_RNA.y.y.y	percent.mt.y	nCount_Protein.y	nFeature_Protein.y	T_clonotype_id.y	T_cdr3s_aa.y	B_clonotype_id.y	B_cdr3s_aa.y	scrublet_score.y	integrated_snn_res.0.1.y	integrated_snn_res.0.15.y	seurat_clusters.y	integrated_snn_res.0.12.y	percent.ribro.y	percent.Hb.y	S.Score.y	G2M.Score.y	Phase.y	old.ident.y	patient.y	TYPE.y	sample.y	sex.y	batch.y	age.y	condition.y	njoint.y	ana.y	global1.y	onset.y	maxjoint0623.y	label15.y	simple15.y	clusters2312.y	oldclusters.y	label15s.y	subclusters.y	main_clusters.y	percent.mito.y	humap_fgraph_res.1.y	humap_fgraph_res.0.8.y	cell.states.y	cell.states.knn.prob.y	donorID.y	tissue.y	cell_type_pred_knn	cell_type_pred_knn_prob
	<chr>	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>	<int>	<chr>	<int>	<int>	<dbl>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<int>	<int>	<fct>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<chr>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<int>	<chr>	<chr>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<fct>	<fct>	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>	<int>	<chr>	<int>	<int>	<dbl>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<int>	<int>	<fct>	<int>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<chr>	<int>	<chr>	<chr>	<chr>	<chr>	<dbl>	<chr>	<int>	<chr>	<chr>	<int>	<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<fct>	<fct>	<fct>	<dbl>	<int>	<chr>	<fct>	<dbl>
1	AAAGATGAGCAGCGTA-1_1	JIA_Croft	1350	760	AAAGATGAGCAGCGTA-1_1	JIA_Croft	1350	760	PBMC1	1350	760	3.1111111	340	99	clonotype4891	TRB:CASSLAIRNYGYTF	NA	NA	0.06086957	1	1	0	1	0.1507320	0	0.01850598	0.068792887	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells	3.1111111	8	0	Naive Treg	0.4333333	91	Blood	JIA_Croft	1350	760	AAAGATGAGCAGCGTA-1_1	JIA_Croft	1350	760	PBMC1	1350	760	3.1111111	340	99	clonotype4891	TRB:CASSLAIRNYGYTF	NA	NA	0.06086957	1	1	0	1	0.1507320	0	0.01850598	0.068792887	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ KLRB1+ memory T	CD4+ KLRB1+ memory	CD4+ T cells	CD4+ KLRB1+ memory T	CD4+ T cells	3.1111111	8	0	Naive Treg	0.4333333	91	Blood	Naive Treg	0.4333333
2	AAAGATGTCCTAGTGA-1_1	JIA_Croft	1485	773	AAAGATGTCCTAGTGA-1_1	JIA_Croft	1485	773	PBMC1	1485	773	1.4814815	292	93	clonotype77	TRA:CAVRLSWGKLQF	NA	NA	0.08163265	1	1	0	1	0.1495413	0	-0.04059452	0.008258175	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.4814815	8	0	Naive Treg	0.5666667	91	Blood	JIA_Croft	1485	773	AAAGATGTCCTAGTGA-1_1	JIA_Croft	1485	773	PBMC1	1485	773	1.4814815	292	93	clonotype77	TRA:CAVRLSWGKLQF	NA	NA	0.08163265	1	1	0	1	0.1495413	0	-0.04059452	0.008258175	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.4814815	8	0	Naive Treg	0.5666667	91	Blood	CD25int Treg	0.5666667
3	AAAGTAGCACTAGTAC-1_1	JIA_Croft	1567	847	AAAGTAGCACTAGTAC-1_1	JIA_Croft	1567	847	PBMC1	1567	847	2.1059349	330	100	clonotype1715	TRB:CASSSQNRAEQYF	NA	NA	0.05263158	1	1	1	1	0.1298440	0	0.01361041	0.082733395	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	2.1059349	8	1	Naive Treg	0.6000000	91	Blood	JIA_Croft	1567	847	AAAGTAGCACTAGTAC-1_1	JIA_Croft	1567	847	PBMC1	1567	847	2.1059349	330	100	clonotype1715	TRB:CASSSQNRAEQYF	NA	NA	0.05263158	1	1	1	1	0.1298440	0	0.01361041	0.082733395	G2M	G2M	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	2.1059349	8	1	Naive Treg	0.6000000	91	Blood	Naive Treg	0.6333333
4	AACCATGGTCGCATCG-1_1	JIA_Croft	2589	1048	AACCATGGTCGCATCG-1_1	JIA_Croft	2589	1048	PBMC1	2589	1048	0.9269988	278	89	clonotype751	TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF	NA	NA	0.15625000	1	1	1	1	0.1489572	0	0.04661555	-0.045203442	S	S	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	0.9269988	8	1	Naive Treg	0.9333333	91	Blood	JIA_Croft	2589	1048	AACCATGGTCGCATCG-1_1	JIA_Croft	2589	1048	PBMC1	2589	1048	0.9269988	278	89	clonotype751	TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF	NA	NA	0.15625000	1	1	1	1	0.1489572	0	0.04661555	-0.045203442	S	S	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	0.9269988	8	1	Naive Treg	0.9333333	91	Blood	Naive Treg	0.9333333
5	AACCATGGTTGCTCCT-1_1	JIA_Croft	2049	1102	AACCATGGTTGCTCCT-1_1	JIA_Croft	2049	1102	PBMC1	2049	1102	1.2201074	279	90	clonotype9315	TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF	NA	NA	0.03914591	1	1	6	1	0.1104492	0	-0.04041617	-0.052987582	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.2201074	8	6	CD25int Treg	0.5666667	91	Blood	JIA_Croft	2049	1102	AACCATGGTTGCTCCT-1_1	JIA_Croft	2049	1102	PBMC1	2049	1102	1.2201074	279	90	clonotype9315	TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF	NA	NA	0.03914591	1	1	6	1	0.1104492	0	-0.04041617	-0.052987582	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CD4+ FOXP3+ Tregs	CD4+ FOXP3+ Tregs	CD4+ T cells	CD4+ FOXP3+ Tregs	CD4+ T cells	1.2201074	8	6	CD25int Treg	0.5666667	91	Blood	ISG Treg	0.5333333
6	AACGTTGAGGATGGAA-1_1	JIA_Croft	4757	1769	AACGTTGAGGATGGAA-1_1	JIA_Croft	4757	1769	PBMC1	4757	1769	1.8499054	293	92	clonotype1286	TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF	NA	NA	0.09497207	1	1	1	1	0.1185773	0	-0.02773367	-0.003626403	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	1.8499054	8	1	Naive Treg	0.7666667	91	Blood	JIA_Croft	4757	1769	AACGTTGAGGATGGAA-1_1	JIA_Croft	4757	1769	PBMC1	4757	1769	1.8499054	293	92	clonotype1286	TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF	NA	NA	0.09497207	1	1	1	1	0.1185773	0	-0.02773367	-0.003626403	G1	G1	91	Blood	B91	F	batch1	6.9	P-Oligo	1	N	T cell	518	1	CD4+ T	T cells	CXCL13+ TpH	CXCL13+ TpH	CD4+ T cells	CXCL13+ TpH	CD4+ T cells	1.8499054	8	1	Naive Treg	0.7666667	91	Blood	Naive Treg	0.7333333

## Add umap and cell state predictions to seurat obj.
tregs@reductions[["umap"]] <- CreateDimReducObject(embeddings = query$umap[rownames(tregs@meta.data),], key = "UMAP_", assay = DefaultAssay(tregs))
tregs <- AddMetaData(tregs, mt$cell_type_pred_knn, "cell.states")
tregs <- AddMetaData(tregs, mt$cell_type_pred_knn_prob, "cell.states.knn.prob")

Save Symphony obj

# saveRDS(query, "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony.RDS")
# saveRDS(tregs, "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony_Seurat.RDS")


tregs <- readRDS("/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony_Seurat.RDS")
query <- readRDS( "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony.RDS")

tregs

An object of class Seurat 
30803 features across 5523 samples within 1 assay 
Active assay: RNA (30803 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 4 dimensional reductions calculated: pca, harmony, humap, umap

Further exploration

query$meta_data %>% names()

1. 'orig.ident'
2. 'nCount_RNA'
3. 'nFeature_RNA'
4. 'cellID'
5. 'orig.ident.x'
6. 'nCount_RNA.x'
7. 'nFeature_RNA.x'
8. 'orig.ident.y'
9. 'nCount_RNA.y'
10. 'nFeature_RNA.y'
11. 'percent.mt'
12. 'nCount_Protein'
13. 'nFeature_Protein'
14. 'T_clonotype_id'
15. 'T_cdr3s_aa'
16. 'B_clonotype_id'
17. 'B_cdr3s_aa'
18. 'scrublet_score'
19. 'integrated_snn_res.0.1'
20. 'integrated_snn_res.0.15'
21. 'seurat_clusters'
22. 'integrated_snn_res.0.12'
23. 'percent.ribro'
24. 'percent.Hb'
25. 'S.Score'
26. 'G2M.Score'
27. 'Phase'
28. 'old.ident'
29. 'patient'
30. 'TYPE'
31. 'sample'
32. 'sex'
33. 'batch'
34. 'age'
35. 'condition'
36. 'njoint'
37. 'ana'
38. 'global1'
39. 'onset'
40. 'maxjoint0623'
41. 'label15'
42. 'simple15'
43. 'clusters2312'
44. 'oldclusters'
45. 'label15s'
46. 'subclusters'
47. 'main_clusters'
48. 'percent.mito'
49. 'humap_fgraph_res.1'
50. 'humap_fgraph_res.0.8'
51. 'cell.states'
52. 'cell.states.knn.prob'
53. 'donorID'
54. 'tissue'
55. 'cell_type_pred_knn'
56. 'cell_type_pred_knn_prob'

fig.size(6,6)
count.groups <- query$meta_data %>% group_by(patient) %>% mutate(totalP = n()) %>% ungroup()%>% 
select(cell_type_pred_knn, TYPE, patient, totalP) %>% 
group_by(cell_type_pred_knn, TYPE, patient) %>% mutate(nType = n()) %>%
mutate(propType = nType/ totalP) %>% ungroup() %>% 
group_by(cell_type_pred_knn, TYPE) %>% summarize(prop.average = mean(propType))
count.groups
ggplot(count.groups, aes(x = TYPE, y = prop.average, fill = cell_type_pred_knn)) +
       geom_col(position = "fill") + scale_fill_manual(values = cell.state.colors.v2$cell.states) + theme_base(base_size = 20)

`summarise()` has grouped output by 'cell_type_pred_knn'. You can override using the `.groups`
argument.

A grouped_df: 29 × 3

cell_type_pred_knn	TYPE	prop.average
<fct>	<chr>	<dbl>
Naive Treg	Blood	0.098614599
Naive Treg	SF	0.009955783
Naive Treg	Tissue	0.016194000
CD25int Treg	Blood	0.084877277
CD25int Treg	SF	0.098925497
CD25int Treg	Tissue	0.112756653
CD25high Treg	Blood	0.049476911
CD25high Treg	SF	0.290763208
CD25high Treg	Tissue	0.200442745
CD25highCXCR6+ Treg	Blood	0.002112659
CD25highCXCR6+ Treg	SF	0.191614706
CD25highCXCR6+ Treg	Tissue	0.158955875
AREG+ Treg	Blood	0.005190914
AREG+ Treg	SF	0.023143855
AREG+ Treg	Tissue	0.205168262
TNFAIP3+ Treg	SF	0.006382039
TNFAIP3+ Treg	Tissue	0.135917068
CD161+mem. Treg	Blood	0.008174522
CD161+mem. Treg	SF	0.015064627
CD161+mem. Treg	Tissue	0.017886946
ISG Treg	Blood	0.006153404
ISG Treg	SF	0.070800637
ISG Treg	Tissue	0.007396221
GZM Treg	Blood	0.004656989
GZM Treg	SF	0.004756350
GZM Treg	Tissue	0.007300054
Prolif.	Blood	0.001792115
Prolif.	SF	0.002680202
Prolif.	Tissue	0.005463594

fig.size(6,6)
count.groups <- query$meta_data %>% select(cell_type_pred_knn, condition) %>% group_by(cell_type_pred_knn, condition) %>% count()


fig.size(6,6)
count.groups <- query$meta_data %>% group_by(patient) %>% mutate(totalP = n()) %>% ungroup()%>% 
select(cell_type_pred_knn, condition, patient, totalP) %>% 
group_by(cell_type_pred_knn, condition, patient) %>% mutate(nType = n()) %>%
mutate(propType = nType/ totalP) %>% ungroup() %>% 
group_by(cell_type_pred_knn, condition) %>% summarize(prop.average = mean(propType))
 
ggplot(count.groups, aes(x = condition, y = prop.average, fill = cell_type_pred_knn)) +
    geom_col(position = "fill") + scale_fill_manual(values = cell.state.colors.v2$cell.states) + 
    theme_base(base_size = 20) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(count.groups, aes(x = condition, y = prop.average, fill = cell_type_pred_knn)) +
    geom_col(position = "stack") + scale_fill_manual(values = cell.state.colors.v2$cell.states) + 
    theme_base(base_size = 20) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

`summarise()` has grouped output by 'cell_type_pred_knn'. You can override using the `.groups`
argument.

table(query$meta_data$condition, query$meta_data$patient)

          
            84  88  91  95  96  97  98  99 811 814 817 910
  Ex-Oligo   0   0   0   0   0   0   0   0   0   0   0 552
  P-Oligo    0   0 558 592   0   0 820 270 281   0 496   0
  Poly-neg 563   0   0   0   0 263   0   0   0   0   0   0
  PsA        0   0   0   0 194   0   0   0   0 238   0   0
  RF+Poly    0 696   0   0   0   0   0   0   0   0   0   0

table(query$meta_data$cell_type_pred_knn, query$meta_data$patient)

                     
                       84  88  91  95  96  97  98  99 811 814 817 910
  Naive Treg           50  38  73  90   4  42  54   2   9  33  16  11
  CD25int Treg        139  88 100 146  22  29  68  22  31  48  33  98
  CD25high Treg       183 207 200 171  96  59 337  74 185  44 168  87
  CD25highCXCR6+ Treg  23 195 102  96  40  85 202  57  36  45 156  33
  AREG+ Treg           84  74   5  56  16   7  24  61   2  24  31 195
  TNFAIP3+ Treg        41  66   1   7   4   3  68  47   2  34  64 118
  CD161+mem. Treg       4  11  12  11   8   8  18   3   9   5  16   5
  ISG Treg             37  16  60   9   2  30  45   1   6   5   7   1
  GZM Treg              1   0   4   6   0   0   2   2   1   0   3   4
  Prolif.               1   1   1   0   2   0   2   1   0   0   2   0

obj$condition %>% unique()

1. 'P-Oligo'
2. 'Poly-neg'
3. 'RF+Poly'
4. 'PsA'
5. 'Ex-Oligo'

fig.size(6,6)
query$meta_data <- query$meta_data %>% group_by(patient, T_clonotype_id) %>% mutate(cloneSize = n()) %>% 
mutate(cloneSize = if_else(is.na(T_clonotype_id), 0, cloneSize))%>% 
ungroup()

qp <- query$umap %>% as.data.frame %>% cbind(query$meta_data) %>% arrange(cloneSize) %>% 
ggplot(aes(x=UMAP1, y=UMAP2, color=as.factor(cloneSize))) + geom_point(size = 0.1, position = ) + 
theme_bw(base_size = 18) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +coord_fixed()


qp

metadata$cloneSize %>% table

.
   0    1    2    3    4    5    6    7    8    9   12   23   26 
 340 3603  822  297  136  100   48   49   40   27   12   23   26 

fig.size(6,8)
metadata <- query$meta_data %>% 
mutate(cloneCategory = ifelse(is.na(T_clonotype_id), "No TCR information",
                                ifelse(cloneSize == "1", "Singlet",
                                ifelse(cloneSize %in% c(2:3), "Small (1 < X <= 3)",
                                       ifelse(cloneSize %in% c(4:10), "Medium (3 < X <= 10)",
                                            ifelse(cloneSize %in% c(11:30), "Large (10 < X <= 30)", 
                                                     "Hyperexpanded (30 < X)")))))) %>%
    mutate(clonal = ifelse(cloneSize > 1, "Clonal", "Not Clonal"))
# metadata$cloneCategory <- tidyr::replace_na(metadata$cloneCategory, "No TCR information")


metadata$cloneCategory <- factor(metadata$cloneCategory, levels = c("Hyperexpanded (30 < X)", 
                                                                "Large (10 < X <= 30)", 
                                                                "Medium (3 < X <= 10)", 
                                                                "Small (1 < X <= 3)", 
                                                                "Singlet", 
                                                                "No TCR information"))

category.order <- c("Hyperexpanded (30 < X)", "Large (10 < X <= 30)",
                  "Medium (3 < X <= 10)","Small (1 < X <= 3)",
                  "Singlet", "No TCR information", "No TCR data")

qp <- query$umap %>% as.data.frame %>% cbind(metadata) %>%  
ggplot(aes(x=UMAP1, y=UMAP2, color=cloneCategory)) + geom_point(size = 1) + 
theme_bw(base_size = 18) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +coord_fixed() +
scale_color_manual(values = colors.no.info)


qp

fig.size(5,8)
metadata %>%
    group_by(cell_type_pred_knn, cloneCategory) %>%
    dplyr::count() %>%
    group_by(cell_type_pred_knn) %>%
    mutate(Percent=100*n/sum(n)) %>%
    ungroup() %>%
    ggplot(aes(x=cell_type_pred_knn, y=Percent, fill=cloneCategory)) +
        geom_col() +
        # scale_fill_brewer(palette = "Blues", direction = -1, name = "Clone Type") +
        scale_fill_manual(values = colors.no.info) +
        scale_y_continuous(labels = function(x) paste0(x, "%")) +
        ggtitle(NULL) +
        labs(fill = "Sample") +
        xlab("cell states") + ylab(NULL) +
        theme_classic() + 
        theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
              legend.text = element_text(size=12),
              legend.title = element_text(size = 14),
            plot.margin = margin(t = 0, r = 0, b = 0, l = 2, unit = "cm"))  


metadata %>%
    group_by(TYPE, cloneCategory) %>%
    dplyr::count() %>%
    group_by(TYPE) %>%
    mutate(Percent=100*n/sum(n)) %>%
    ungroup() %>%
    ggplot(aes(x=TYPE, y=Percent, fill=cloneCategory)) +
        geom_col() +
        # scale_fill_brewer(palette = "Blues", direction = -1, name = "Clone Type") +
        scale_fill_manual(values = colors.no.info) +
        scale_y_continuous(labels = function(x) paste0(x, "%")) +
        ggtitle(NULL) +
        labs(fill = "Sample") +
        xlab("cell states") + ylab(NULL) +
        theme_classic() + 
        theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
              legend.text = element_text(size=12),
              legend.title = element_text(size = 14),
            plot.margin = margin(t = 0, r = 0, b = 0, l = 2, unit = "cm"))  
# merged@meta.data %>%
#     group_by(tissue, cloneCategory) %>%
#     dplyr::count() %>%
#     group_by(tissue) %>%
#     mutate(Percent=100*n/sum(n)) %>%
#     ungroup() %>%
#     ggplot(aes(x=tissue,y=Percent, fill=cloneCategory)) +
#         geom_col() +
#         # scale_fill_brewer(palette = "Blues", direction = -1, name = "Clone Type") +
#         scale_fill_manual(values = colors.no.info) +
#         scale_y_continuous(labels = function(x) paste0(x, "%")) +
#         ggtitle(NULL) +
#         labs(fill = "Sample") +
#         xlab("tissue") + ylab(NULL) +
#         theme_classic() + 
#         theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
#               legend.text = element_text(size=12),
#               legend.title = element_text(size = 14))

fig.size(8, 12)
#getting clones for each cluster
# unique.states <- unique(merged$cell.states)
unique.tissue <- unique(metadata$tissue)
all.states.clones <- list()
for (s in unique.tissue){
        a <- metadata %>% filter(tissue == s) %>% 
    select(T_clonotype_id) %>% na.omit() %>% distinct(T_clonotype_id)
        all.states.clones[s] = a
    }

#using the upsetR package
# all.states.clones 
UpSetR::upset(UpSetR::fromList(all.states.clones), sets = unique.tissue, keep.order = F, 
              mainbar.y.label = "Number of Overlapping Clones", text.scale = 2,
              sets.x.label = "Unique Clones in tissue", order.by = "freq")

fig.size(8, 18)
#getting clones for each cluster
# unique.states <- unique(merged$cell.states)
unique.tissue <- unique(metadata$cell_type_pred_knn)
all.states.clones <- list()
for (s in unique.tissue){
        a <- metadata %>% filter(cell_type_pred_knn == s) %>% 
    select(T_clonotype_id) %>% na.omit() %>% distinct(T_clonotype_id)
        all.states.clones[s] = a
    }

#using the upsetR package
# all.states.clones 
UpSetR::upset(UpSetR::fromList(all.states.clones), sets = unique.tissue, keep.order = F, 
              mainbar.y.label = "Number of Overlapping Clones", text.scale = 2,
              sets.x.label = "Unique Clones in cell.state", order.by = "freq")

Markers expression

1

fig.size(6,10)
Idents(tregs) <- "cell.states"
DimPlot(tregs, reduction = "umap", cols = cell.state.colors.v2$cell.states, pt.size = 0.8)

fig.size(8, 12)
query$umap %>% as.data.frame %>% 
cbind(query$meta_data) %>% 
ggplot(aes(x=UMAP1, y=UMAP2, color=cell.states)) + geom_point() + 
scale_color_manual(values = cell.state.colors.v2$cell.states) +
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) + 
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + coord_fixed()

fig.size(3,4)
FeaturePlot(tregs, features = "IL2RA", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "FOXP3", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "AREG", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "CXCR6", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "CXCR4", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "CXCL13", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "TCF7", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "NR3C1", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "NR4A2", reduction = "umap") + coord_fixed()

#non expressed - EGF, EPGN, NRG3, 
ligands <- c("VEGF", "AREG", "TGFA", "HBEGF", "EREG", "BTC", "NRG4", "NRG1", "NRG2")
ligands2 <- c("NR3C1", "NR4A2", "NR4A1")

fig.size(10,10)
FeaturePlot(tregs, features = ligands, reduction = "umap") 

fig.size(3,10)
FeaturePlot(tregs, features = ligands2, reduction = "umap", ncol =3)

Warning message:
“The following requested variables were not found: VEGF”

fig.size(5,15)
FeaturePlot(tregs, features = "AREG", reduction = "umap", split.by = "TYPE")

Tregs.markers <- wilcoxauc(tregs, "cell.states")

options(repr.matrix.max.cols=150, repr.matrix.max.rows=200)

Tregs.markers %>% 
    group_by(group) %>% 
    filter(padj < 0.05 & auc > 0.6) %>% 
    arrange(group, desc(logFC)) %>% 
    #slice_max(n = 40, order_by = logFC)%>% 
    dplyr::select(feature, group) %>% 
    mutate(row = row_number()) %>% 
    tidyr::pivot_wider(names_from = group, values_from = feature) -> Treg.show
Treg.show[1:50,]

A tibble: 50 × 11

row	AREG+ Treg	CD161+mem. Treg	CD25high Treg	CD25highCXCR6+ Treg	CD25int Treg	GZM Treg	ISG Treg	Naive Treg	Prolif.	TNFAIP3+ Treg
<int>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
1	FOS	AQP3	LGALS1	TNFRSF18	RPS29	GZMK	ISG15	LEF1	DEK	CREM
2	JUNB	GIMAP7	S100A4	TNFRSF4	IFITM1	CCL5	MX1	EEF1B2	STMN1	JUNB
3	DUSP1	GPR25	LGALS3	CTSC	RPS12	GZMA	IFI6	TCF7	TUBA1B	ZFP36
4	ZNF331	CD52	S100A6	LINC01943	LEF1	IFNG-AS1	LY6E	RPS3A	BCAP31	FOSB
5	JUND	S100A11	S100A10	CTLA4	RPS27	GIMAP5	IFIT1	RPS12	NME2	DUSP2
6	FOSB	KLRB1	ANXA2	BATF	TXNIP	DENR	ISG20	RPL3	DUT	RGCC
7	JUN	CNN2	EMP3	CD74	RPL32	PRKCQ-AS1	HERC5	RPS6	NDUFA6	DUSP1
8	CREM	FOXP1	CD74	TIGIT	RPL30	GLMN	EPSTI1	TXNIP	TXN	ZNF331
9	ZFP36	GIMAP4	HLA-DRB1	SPOCK2	RPL13	TIMM22	XAF1	RPL13	CCDC167	FOS
10	DUSP2	NOSIP	S100A11	PKM	RPS8	CHI3L2	STAT1	RPS5	YBX1	TNFAIP3
11	NR4A2	GIMAP5	CRIP1	CXCR6	RPS18	ACACA	OAS1	RPS8	GAPDH	CXCR4
12	CD69	TRAF3IP3	ISG20	IL2RA	RPS21	SREBF1	IFITM1	RPL32	SUMO3	DNAJB1
13	TNFAIP3	EMP3	HLA-DPA1	CD2	RPS3A	KRBOX4	MX2	RPS13	PFN1	LMNA
14	RGCC	TLE5	SELPLG	TYMP	RPL34	NA	SP100	SELL	C17orf49	NR4A2
15	DNAJB1	LIMS1	VIM	GAPDH	RPS13	NA	IFI16	CCR7	PRDX3	JUND
16	CXCR4	SLC25A5	HLA-DPB1	ENO1	RPS28	NA	OASL	RPL10A	GSTP1	JUN
17	SLC2A3	RPLP0	LY6E	FOXP3	RPL8	NA	SMCHD1	RPS18	TWF2	PHLDA1
18	KLF6	COTL1	CD52	TNFRSF1B	RPL10	NA	IRF7	RPLP2	DCTPP1	DNAJA1
19	DNAJA1	KLF3	CLIC1	HLA-DRB1	RPL3	NA	TRIM22	EEF1G	CDT1	CD69
20	BTG2	ARPC1B	CD99	LGALS3	RPL18A	NA	EIF2AK2	RPL5	PHF19	PPP1R15A
21	YPEL5	LAPTM5	SYNE2	SRGN	RPL11	NA	RNF213	RPL9	OVCA2	DUSP4
22	PPP1R15A	ACTG1	SH3BGRL3	CARD16	RPL39	NA	SAMD9	GIMAP7	ACTB	YPEL5
23	RGS1	MGAT4A	LIME1	PGAM1	RPL36	NA	CMPK2	RPL29	PSMA3	TSC22D3
24	SOCS1	RCSD1	IL32	RNF213	RPS3	NA	IFI44L	FOXP1	MDH1	TENT5C
25	TENT5C	RPL7	TLE5	LGALS1	RPLP2	NA	UBE2L6	RPS23	TRBV6-1	SELENOK
26	ZFP36L2	IKZF3	LSP1	TPI1	RPL22	NA	SLFN5	RPL21	HLA-DMA	HSPH1
27	TSC22D3	RPSA	CYBA	IL2RB	RPL21	NA	DRAP1	RPL4	SNRPE	KLF6
28	FTH1	MSN	TMSB10	IL2RG	RPL14	NA	LGALS9	RPL7	RBPJ	BTG2
29	UBE2S	RIPOR2	HLA-A	CST7	RPL19	NA	SP110	RPL22	MCM5	SLC2A3
30	SELENOK	COX6C	HLA-C	MYL6	RPS23	NA	SAMD9L	RPL18A	PSMB2	RGS1
31	NR3C1	CCR4	PFN1	LAG3	RPS6	NA	TAP1	RPS21	MCM7	FTH1
32	HSPA1B	ABRACL	ARHGDIB	DUSP4	RPS27A	NA	PARP9	RPL13A	PDAP1	NR4A1
33	RGS2	S1PR4	TMSB4X	PHACTR2	EEF1B2	NA	BST2	EEF1A1	COA3	SRGN
34	HSP90AA1	CYTH1	B2M	TIMP1	RPL41	NA	IFIT3	RPL30	SUCLG1	SAMSN1
35	HSPH1	RORA	NA	CKLF	RPL12	NA	OAS3	RPL36A	PPIA	HSP90AA1
36	TAGAP	CCND3	NA	IGFLR1	RPL5	NA	USP18	RPL34	MRPL22	UBE2S
37	SRGN	RPL4	NA	PRDM1	RPL10A	NA	SPATS2L	RPS27A	DDX49	SOCS1
38	TSPYL2	CDC25B	NA	LAYN	RPL37	NA	IFI35	RPL38	SLBP	MYADM
39	FAAH2	CD6	NA	PTPRJ	RPS15A	NA	IFI44	RPS2	CISD1	BTG1
40	CCNH	MKRN1	NA	S100A4	RPL29	NA	OAS2	GAS5	RPA3	CSRNP1
41	H3F3B	MVP	NA	DNPH1	FAU	NA	SYNE2	RPL39	CLSPN	NR4A3
42	NR4A1	SYTL1	NA	ACP5	RPS14	NA	ADAR	RPS3	PTTG1	NAMPT
43	PHLDA1	PLSCR3	NA	IL12RB2	RPS5	NA	TYMP	RPS14	SLC35E4	HSPA1B
44	GEM	ITGB7	NA	CD4	RPL38	NA	RSAD2	RPL19	CMC1	CCNH
45	TUBB4B	GIMAP1	NA	CBLB	RPS25	NA	ZBP1	RPL36	POLR3K	PTPN22
46	LMNA	HMOX2	NA	CACYBP	FTL	NA	PSMB8	LDHB	PCLAF	FOSL2
47	MTRNR2L12	AC004687.1	NA	SH2D2A	RPL35A	NA	USP15	RPL37	APOBEC3G	HSPA8
48	DDIT4	TC2N	NA	TXN	RPL26	NA	GBP1	RPS20	NABP2	IRF1
49	EIF1	CDC123	NA	GBP5	RPL18	NA	FOXP3	RPL23A	LIG1	NFKBIA
50	SOCS3	TMSB10	NA	UCP2	RPS7	NA	DDX60L	PABPC1	MRPS7	HSP90AB1

Metabolism Genes

mg <- read.csv("/data/brennerlab/Shani/projects/Treg/useful_tables/metabolism_genes_manually_curated.csv")
mg %>% head

A data.frame: 6 × 3

	Glycolysis	FAO	OXPHOS
	<chr>	<chr>	<chr>
1	HK1	ACADVL	PDHA1
2	HK2	ACADM	DLST
3	PFKP	HADHA	SDHA
4	PGAM1	HADHB	CS
5	ENO1	CPT1A	ACO1
6	PKM	CPT2	IDH1/2/3

fig.size(4,8)
FeaturePlot(tregs, features = mg$Glycolysis, reduction = "umap") + coord_fixed()

Warning message:
“The following requested variables were not found:  HK2,  PFKP,  PGAM1,  ENO1,  PKM,  LDHA,  PDHA1,  TPI1,  GAPDH”

fig.size(8, 12)
FeaturePlot(tregs, features = mg$FAO, reduction = "umap")

Warning message:
“The following requested variables were not found: ”

fig.size(6, 10)
FeaturePlot(tregs, features = mg$OXPHOS, reduction = "umap")

Warning message:
“The following requested variables were not found: ACO1 , IDH1/2/3, OGDH , SUCLA2/SUCLG1/SUCLG2, FH , MDH1/2, DLD ”

Try harmonizid CD and RA Tregs

RA <- readRDS("/data/brennerlab/Shani/projects/Treg/analysis/integrated/integrated.Tregs.rds")

length(rownames(RA))
length(rownames(tregs))

inters <- intersect(rownames(RA), rownames(tregs))
length(inters)

38224

30803

30803

tregs@meta.data <- tregs@meta.data %>% mutate(donorID = patient) %>% mutate(tissue = TYPE)

merged <- merge(RA[inters,], tregs[inters,])

merged

An object of class Seurat 
30803 features across 24432 samples within 1 assay 
Active assay: RNA (30803 features, 0 variable features)
 2 layers present: counts, data

merged <- NormalizeData(merged, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>% 
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>% 
    ScaleData() %>% 
    RunPCA(verbose = F)

Centering and scaling data matrix

unique(merged$orig.ident)
unique(merged$tissue)

1. 'AMP2'
2. 'AMPrep'
3. 'SF.BL'
4. 'HD_Luo'
5. 'JIA_Croft'

1. 'RA.Syn.Tissue'
2. 'RA.PBMC'
3. 'RA.Syn.Fluid'
4. 'HD.PBMC'
5. 'Blood'
6. 'SF'
7. 'Tissue'

fig.size(4,4)
set.seed(1)
merged <- RunHarmony(merged, group.by.vars=c("orig.ident", "donorID", "tissue"), reduction.use = "pca", 
                 plot_convergence = TRUE, max_iter = 10,
                 early_stop = T)

Transposing data matrix

Initializing state using k-means centroids initialization

Harmony 1/10

Harmony 2/10

Harmony 3/10

Harmony converged after 3 iterations

set.seed(1)
merged <- Run_uwot_umap(merged)
merged <- FindClusters(merged, graph.name = 'humap_fgraph', resolution = 0.8, verbose = TRUE)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 24432
Number of edges: 304539

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.7036
Number of communities: 11
Elapsed time: 2 seconds

merged$tissue <- factor(merged$tissue, levels = c("HD.PBMC", "RA.PBMC", "RA.Syn.Fluid", "RA.Syn.Tissue", 'Blood', 'SF', 'Tissue'))

fig.size(8,8)
DimPlot(merged, label = T, repel = T, pt.size = 0.5, label.size = 6) + 
  scale_color_tableau("Tableau 20") + theme(text = element_text(size = 20))
DimPlot(merged, group.by = "donorID")
DimPlot(merged, group.by = "tissue")
DimPlot(merged, group.by = "orig.ident")

DimPlot(merged, group.by = "cell.states") + scale_color_manual(values = cell.state.colors$cell.states)

fig.size(6, 20)
DimPlot(merged, group.by = "cell.states", split.by = "tissue") + scale_color_manual(values = cell.state.colors$cell.states)

Plot RA Treg reference

fig.size(8, 12)
tregRef$umap$embedding %>% as.data.frame %>% 
cbind(tregRef$meta_data) %>% filter(tissue != "HD.PBMC") %>% 
ggplot(aes(x=UMAP1, y=UMAP2, color=cell.states)) + geom_point() + 
scale_color_manual(values = cell.state.colors$cell.states) +
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) + 
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + coord_fixed()

fig.size(12,16)
tregRef$umap$embedding %>% as.data.frame %>% 
cbind(tregRef$meta_data) %>% 
ggplot(aes(x=UMAP1, y=UMAP2, color=cell.states)) + geom_point() + 
scale_color_manual(values = cell.state.colors$cell.states) + facet_wrap(~tissue) +
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) + 
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + coord_fixed()